#!/usr/bin/perl
#===============================================================================
# foldTimeSeries.pl: script to fold a time series
#===============================================================================

use Math::BigFloat;
use POSIX;

# Command line options
# ==============================================================================
print STDERR "============================================\n";
print STDERR "             foldTimeSeries.pl\n";
print STDERR "============================================\n";
$argc = $#ARGV+1;
if ($argc != 4) {
    print STDERR "ERROR: incorrect number of command line args ($argc/4)\n";
    print STDERR "ERROR:\n";
    print STDERR "USAGE: foldTimeSeries.pl [data file] [fold period] [n bins] [output file root name]\n";
    print STDERR "\n";
    exit(99);
}

# Save command line arguments and initialise/declare variables
# ==============================================================================
$dataFile   = $ARGV[0]; 
$foldPeriod = Math::BigFloat->new($ARGV[1]);
$outputFile = $ARGV[3];
$nBins      = $ARGV[2]; # Number of phase bins

$foldFreq   = 1.0/$foldPeriod;
$binWidth   = 1.0/$nBins; # Bin width in phase units (0->2pi)

$epoch      = 2450023.62;                           # reference epoch
$epoch      = Math::BigFloat->new(0.0);             # no correction
$epoch      = 2450023.62 - 500.0 * $foldPeriod;     # speeds2004 data
$epoch      = 2450023.62 - 650.0 * $foldPeriod;     #
$epoch      = 2450023.62 - 3300.0 * $foldPeriod;    #  
$epoch      = 2450023.62 - 485.0 * $foldPeriod;     # speeds data 

$outputFile = "${outputFile}-p${foldPeriod}-n${nBins}.dat";



# Read the data file
# ==============================================================================
$i = 0;
$nData = 0;
$phase      = Math::BigFloat->new(0.0);      # Phase of the signal (in units of 2pi)
open(DATA,"$dataFile") || die "cannot open: $dataFile\n\n";
while(<DATA>) {
    split;
    $time[$i] = Math::BigFloat->new($_[0]);
    $amp[$i] = Math::BigFloat->new($_[1]);
    $a = Math::BigFloat->new(0.0);
    $b = Math::BigFloat->new(0.0);
    $a =  ($time[$i]-$epoch)/$foldPeriod;
    $b =  int(($time[$i]-$epoch)/$foldPeriod);
    $phase[$i] = ($a - $b);
#    printf(STDERR "time:%8.4f\tamp:%8.3f\tphase: %8.3f\n",
#            $time[$i],$amp[$i],$phase[$i]);
    $i++;
}
close(DATA);
$nData = $i;
 
$timeDelta = $time[0] - $epoch;
# Print a summary of the input options
# ==============================================================================
print STDERR "\n";
print STDERR "--------------------------------------------\n";
print STDERR "Datafile        = $dataFile\n";
print STDERR "nData           = $nData\n";
print STDERR "epoch           = $epoch\n";
print STDERR "start time      = $time[0]\n";
printf(STDERR "time delta      = %.6f\n", $timeDelta);
printf(STDERR "Folding freq    = %.6f\n",$foldFreq);
print STDERR "Folding period  = $foldPeriod\n";
print STDERR "nBins           = $nBins\n";
print STDERR "Bin width       = $binWidth (x 2pi)\n";
print STDERR "Output file     = $outputFile\n";
print STDERR "--------------------------------------------\n";
print STDERR "\n";

# Put data into bins
# ==============================================================================
$totalBinned = 0;
$binOffset = Math::BigFloat->new($binWidth/2.0);
for ($iBin = 0;$iBin < $nBins; $iBin++) {
    $binStart[$iBin] = Math::BigFloat->new(0.0);
    $binCentre[$iBin]= Math::BigFloat->new(0.0);
    $binEnd[$iBin]   = Math::BigFloat->new(0.0);
    $binStart[$iBin] = $iBin * $binWidth;
    $binCentre[$iBin]= $binStart[$iBin] + $binOffset;
    $binEnd[$iBin]   = $binStart[$iBin] + $binWidth;
    $binCount[$iBin] = 0;
    $binTotal[$iBin] = Math::BigFloat->new(0.0);
    $binMean[$iBin]  = Math::BigFloat->new(0.0);
#    printf(STDERR "[bin %2d]\t($binStart[$iBin] <= x < $binEnd[$iBin])\n",$iBin);
    for ($i=0;$i<$nData;$i++) {
        if ( ($phase[$i] >= $binStart[$iBin]) && ($phase[$i] < $binEnd[$iBin]) ) {
            $binCount[$iBin] += 1;
            $binTotal[$iBin] += $amp[$i];
            $totalBinned += 1;
#            printf(STDERR "\t%3d: phase = %.4f (amp = %.3f total = %.3f)\n",$i, $phase[$i], $amp[$i],$binTotal[$iBin]);
        }
    }
    if ($binCount[$iBin] > 0) {
        $binMean[$iBin] = $binTotal[$iBin]/$binCount[$iBin];
    }
#    printf STDERR "[bin %2d]\t(Found $binCount[$iBin] points, mean = %.4f)\n\n",$iBin, $binMean[$iBin]; 
}
if ($totalBinned != $nData) {
    print STDERR "ERROR: not all the data points were binned!\n";
    exit(1);
}

# Caluclate mean amplitude and standard deviation in each bin
# ==============================================================================
for ($iBin=0; $iBin<$nBins; $iBin++) { 
    $binStdDev[$iBin] = Math::BigFloat->new(0.0);
    for ($i = 0; $i < $nData; $i++) {
        if (($phase[$i] >= $binStart[$iBin]) && ($phase[$i] < $binEnd[$iBin]) ) {
            $sumSqDev[$iBin] += ($amp[$i] - $binMean[$iBin])**2;
#            printf(STDERR "amp = %.6f mean %.6f \n", $amp[$i], $binMean[$iBin]);
        }
    }
    if ($binCount[$iBin] > 0) {
        $binStdDev[$iBin] = sqrt($sumSqDev[$iBin]/$binCount[$iBin]);
    }
#    printf(STDERR "($iBin) count = $binCount[$iBin] mean = %.5f sumSq = %.5f stddev = %.5f\n",$binMean[$iBin],$sumSqDev[$iBin],$binStdDev[$iBin]); 
}


open(OUT,">$outputFile") || die "cannot open: $outputFile\n\n";
print STDERR "--------------------------------------------\n";
#========== print output to STDOUT ===========================   
print OUT "# ibin | phase | mean | stdDev | #points \n";
print STDERR "# ibin | phase | mean | stdDev | #points \n";
for ($iBin=0;$iBin<$nBins;$iBin++) { 
    printf(STDOUT "%d %9.3f %9.5f %9.5f %d\n",$iBin, $binCentre[$iBin],$binMean[$iBin],$binStdDev[$iBin],$binCount[$iBin]);
    printf(OUT    "%d %9.3f %9.5f %9.5f %d\n",$iBin, $binCentre[$iBin],$binMean[$iBin],$binStdDev[$iBin],$binCount[$iBin]);
    
}    
print STDERR "--------------------------------------------\n";

close(OUT);
